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FAST AND ACCURATE COMPUTATION OF THE LOGARITHMIC CAPACITY 

OF COMPACT SETS 

JORG LIESEN*, OLIVIER SETEt, AND MOHAMED M.S. NASSER* 


Abstract. We present a numerical method for computing the logarithmic capacity of compact subsets of C, 
which are bounded by Jordan curves and have finitely connected complement. The subsets may have several compo¬ 
nents and need not have any special symmetry. The method relies on the conformal map onto lemniscatic domains 
and, computationally, on the solution of a boundary integral equation with the Neumann kernel. Our numerical 
examples indicate that the method is fast and accurate. We apply it to give an estimate of the logarithmic capacity of 
the Cantor middle third set and generalizations of it. 

Key words, logarithmic capacity, transfinite diameter, Chebyshev constant, conformal map, lemniscatic do¬ 
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AMS subject classifications. 65E05, 30C30, 30C85, 3IA15 

I. Introduction. The logarithmic capacity c{E) of a compact set E in the complex 
plane C is an important invariant that is closely related to polynomial approximation, poten¬ 
tial theory, and conformal mapping. If E is simply connected, then c{E) can be obtained 
via the Riemann map from the complement of E onto the complement of the unit disk. The 
Riemann map is known analytically for some simple sets such as disks, ellipses and intervals. 
This leads to analytic formulas for the logarithmic capacity in terms of the parameters de¬ 
scribing the corresponding sets. Some examples are shown in Table 12.1] below: see also li^ 
p. 135], d p. 557] or d pp. 172-173]. 

The logarithmic capacity can be obtained from the (exterior) Schwarz-Christoffel map 
when E consists of one or several components with polygonal boundaries and has a connected 
complement; see, e.g., |I8] Sections 4.4 and 5.8]. For simply connected sets E this has been 
implemented in the Schwarz-Christoffel Toolbox 111, where the logarithmic capacity is one 
of the outputs of the command extermap. We refer to ii for some interesting examples of 
polygonal sets with special symmetry properties with respect to the real line. 

Considering more complicated sets, in particular the Cantor middle third set, Ransford 
and Rostand pointed out that computing the capacity is “notoriously hard” li^ p. 1499]. They 
derived a method for computing upper and lower bounds for the logarithmic capacity, which 
can in principle be made arbitrarily close to each other. Methods for a direct computation 
of the logarithmic capacity have been proposed in li^ lZTll : also see the survey given in li24l . 
Interest in numerically computing the logarithmic capacity goes back at least to 0. 

In this paper we show that the logarithmic capacity of a fairly wide class of sets can be 
computed fast and accurately using conformal mapping techniques. The method we present 
is based on a conformal map from the complement of E onto a lemniscatic domain, which is 
originally due to Walsh ll3^ and which represents a direct generalization of the Riemann map 
to sets E with several components. We derived a numerical method for computing Walsh’s 
map in ED. The logarithmic capacity is obtained, almost as a by-product, in the first step 
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of that method, and without computing the conformal map itself. Because of the practical 
importance of the logarithmic capacity and the apparent lack of general purpose software for 
its computation, we here derive a stand-alone method and its MATLAB implementation. 

Our method is applicable to any compact set E whose complement in the extended plane 
is connected and bounded by finitely many (sufficiently smooth) Jordan curves. In particular, 
it is required neither that E is connected nor that E has special symmetry properties. Go¬ 
ing beyond the theory presented in ET\ . we place a particular emphasis on the treatment of 
sets E with corners. Moreover, we use a recently developed iterative method for mapping 
parallel slit domains onto domains exterior to ellipses (see Appendix A) in order to com¬ 
pute approximations of the logarithmic capacity of the classical Cantor middle third set and 
generalizations of it. Numerous further numerical examples demonstrate that our method is 
fast and accurate. In our numerical examples with sets for which the logarithmic capacity 
is known analytically, our method typically yields a computed approximation with a relative 
error of order 10“^^. These computations in MATLAB take at most a few seconds on a 
standard laptop. 

The paper is organized as follows. In Section |2] we summarize major facts about the 
logarithmic capacity and its relation to conformal mapping. In Section [3] we describe our 
method for computing the logarithmic capacity and state its MATLAB implementation. In 
Section|4]we give numerical examples. 

2. Background. The transfinite diameter of a compact set E' C C is defined as 

d{E) := lim dn{E) , where dn{E) := max TT \zk — zi\. (2.1) 

n—KXD zi,...,z„e-E 

l<k<£<n 

The existence of the limit in (12.1b was first shown by Fekete ifTOll . who considered the a 
sequence of “generalized diameters” of E. Note that d 2 {E) = l-^i — 22 ] is the 

(usual) diameter of E. In the same article Fekete showed that the transfinite diameter is equal 
to the Chebyshev constant 

t{E) := lim tn, where := min max |z" — p(z)|. (2.2) 

n—¥co deg(p)<n —1 zGE 

Thus, the geometric constant d{E) is closely related to polynomial approximation in the com¬ 
plex plane. This relation can be used to easily show that for both the Chebyshev constant and 
the transfinite diameter it is sufficient to consider compact sets E with a connected comple¬ 
ment. Let E C C be a compact set and let E'^ := C\E denote its complement in the extended 
complex plane C = C U { 00 }. Denote by /C the component of the complement that contains 
the point at infinity (i.e., 00 € /C) and define E = C\/C. Intuitively, E is obtained after 
“filling in the holes” in E. Now the definition of the Chebyshev constant and the maximum 
modulus principle imply that d{E) = t{E) = t(E) = d{E). 

Szego 1321 showed that if the complement of a compact set E is connected and regular 
in the sense that it possesses a Green’s function gE<^ with pole at infinity, then the transfinite 
diameter and the Chebyshev constant are equal to the logarithmic capacitj^ 

c(E) := lim exp(log|z| - 5 ec(z)). (2.3) 

z—¥oo 

Saff 1^ called the result d{E) = t{E) = c{E) \h& fundamental theorem of classical poten¬ 
tial theory. 

'More precisely, Szego showed that d(E) is equal to the Robin constant 7 , which he defined via 
lirriz^oo(log| 2 :| — = log7- In the modem literature the definition of the Robin constant usually gives 

c{E) = exp(- 7 ). 
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E 

c{E) 

disk of radius r 

r 

half-disk of radius r 


ellipse with semi-axes a and b 

\{aEb) 

line segment of length h 


square with side h 

ir(l/4)2/i/7r3/2 

two intervals [—b, —a] U [a, b] 



Table 2.1 

Examples of known logarithmic capacities. 


If = C\E is simply connected (and E is not a single point), then its Green’s function 
is given by gE'=iz) = log|<l>(z)|, where w = is the uniquely determined Riemann 

map from the complement of E to the complement of the unit disk that is normalized by 
$( 00 ) = 00 and $'( 00 ) > 0. Near infinity this map can be written as 

$(z) = ^ + Mo + oQ, 

where l/$'(oo) = p > 0 is called the conformal radius of E. Hence we have 
d{E) = t{E) = c{E) = lim exp(log|z| — log|<i>( 2 ;)|) = g, 

z—^oo 

which shows that the logarithmic capacity can be obtained using (analytical or numerical) 
conformal mapping techniques. The simplest example is a disk E of radius r > 0, for which 
<i>( 2 :) = zjr and thus c(E) = r; see Table 12. ll for some further analytically known examples. 

Walsh proved a direct generalization of the classical Riemann mapping theorem in 
which he replaced the complement of the unit disk by a lemniscatic domain of the form 

i 

£ := {z € C : |(7(2)| > g}, where U{z) := (2-4) 

i=i 

ui,..., € C are pairwise distinct, mi,..., mi are positive real numbers with X]j=i = 

1, and p > 0. A simple calculation shows that log|(7(z)| — log(/r) is the Green’s function 
with pole at infinity for £, so that by ( 12.3b we have 

c(C\/:) = lim exp(log|z| - log|[/( 2 :)| + \og{fi)) = fi. 

z—¥oo 

Walsh proved the following existence theorem in . 

Theorem 2.1. Let E be a compact set whose complement K, = C\E is connected and 
bounded by I Jordan curves. Then there exists a uniquely determined lemniscatic domain C 
of the form dia with p, = c{E) and a uniquely determined conformal map 

^ : K. C with 4>(z) = z + O near infinity. 

The first analytic examples of Walsh’s conformal map onto lemniscatic domains have 
recently been given in . These were applied in OTI in a study of polynomial approxima¬ 
tion problems on disconnected compact sets and, in particular, two real intervals. In II 2 TI we 
developed a numerical method for computing the lemniscatic domain £ and the conformal 
map $ corresponding to a given compact set E with connected complement 1C. 
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3. Computing the logarithmic capacity. As indicated in the formulation of Theo- 
rem l2.ll the logarithmic capacity occurs naturally as one of the parameters defining the lem- 
niscatic domain £ onto which /C = C\E is mapped. Hence the map $ itself is not required 
for computing c{E). A closer inspection of the numerical method developed in 11211 shows 
that the parameters mi,..., m^, /i of £ are indeed computed in a step that can be executed 
separately. Due to the structure of the underlying equations the computation of these param¬ 
eters can be done in a very efficient way. We will now briefly describe this computation. A 
detailed derivation is given in mi. 

As in Theorem 12.11 let E be compact with a finitely connected complement 1C. (As 
noted in Section]!] the connectedness of /C is no restriction for computing the capacity.) We 
assume that the boundary T = dE = dIC of E consists of i Jordan curves Ti,..., T^, which 
satisfy the following smoothness assumption; Each is parameterized by a 27r-periodic 
function rjj : Jj := [0, 27r] —Ej, which is twice continuously differentiable and satisfies 
rij{t) = ^(t) ^ 0 for all t. These assumptions can be relaxed so as to include domains 
with corners, see the precise statement and discussion below. The boundary of E is oriented 
clockwise, so that /C is to the left of the boundary. 

Then a parameterization for the whole boundary E is given by the map 

tCJl, 

; (3.1) 

mit), t e 

where J is the disjoint union of the intervals Ji,..., i.e., J consists of £ copies of [0, 27r]. 

For each j = 1, 2, ...,£, we choose an auxiliary point Uj in the interior of the Jordan 
curve Ej, and define the function 

7j (£) =-log|?7(£) - ajj, i e J. (3.2) 

In practical applications of our method the parameters aj must be specified by the user; cf. 
the MATLAB code shown in Figure [TT] below. Our numerical experience with the method 
suggests that the actual values of the aj are not important, as long as these points are suffi¬ 
ciently far away from the boundary Ej. In the experiments discussed in Section|4|we always 
chose aj close to (or at) the center of the interior of E^. The (blue) dots in Figures |42|and|43] 
show some examples. 

Fet E[ denote the space of all functions / in J, whose restriction to Jj = [0,27r] is a 
real-valued, 27r-periodic and Holder continuous function for each j = 1 , 2 , Define the 
integral operators 

(Nt)(») = 

on H. The kernel of N is called the Neumann kernel. Denoting by I the identity operator on 
El, we can state the following theorem that combines liT\] Theorems 4.1 and 4.3]. 

Theorem 3.1. For each j = 1 , 2 ,... ,£ the integral equation 

(I - N)/i, = -M7, (3.3) 

with as in (13.21) has a unique solution fij € H, and the function 

h , := - (I - N)7,)/2 


(3.4) 
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is real-valued and piecewise constant, that is 


(^) — 5 I' ^ Jk: k — 1 , 2 ,...,-^. 


Furthermore, \og{fi) and the parameters mi,..., mg in (12.41 ) are the unique solution of the 
linear algebraic system 


^1,1 hi2 
^2,1 /i2,2 

hi,e —1 
h2,e — 1 


TOi 

ni2 


1- 

O O . 

1 _ 

hep heg ■ 

hep — 1 


me 


0 

1 1 

1 0 _ 


}og{p)_ 


1 


(3.5) 


This suggests the following method for computing the logarithmic capacity p: 

(1) For each j = 1,... ,i solve the integral equation (13.31) for the unknown function pj. 

(2) Solve the linear algebraic system ( 13.5b of order £ + 1, where the entries hkj in the 
coefficient matrix are computed from ( 13.4b using the known functions yj and the 
functions pj computed in step ([T](- 

The linear algebraic system in step (|2]) is usually quite small and we solve it directly 
using the “backslash” operator in MATLAB. Step ([T]( requires more work: 

As described in 11211 Section 5], the i boundary integral equations ( 13.3b can be solved 
accurately by the Nystrom method with the trapezoidal rule. This method yields a linear al¬ 
gebraic system with a dense nonsymmetric matrix I — B of order In, where n is the number 
of nodes in the discretization of each boundary component. This system can be solved itera¬ 
tively using the GMRES method. Each step of this method requires one multiplication with 
the matrix I — B. Due to the structure of the integral equation, this product can be computed 
efficiently in 0{in) operations using the Fast Multipole Method (FMM). The eigenvalues of 
the matrix I — B are contained in the interval (0, 2] and they cluster around 1. As observed 
in and several other publications (see, e.g., mill), the number of GMRES iterations 
for obtaining a very good approximation of the exact solution is mostly independent of the 
given domain and number of nodes in the discretization of the boundary. 

The method for solving ( 13.3b for the pj and subsequently computing hj in ( 13.4b has been 
implemented in the MATLAB function fbie shown in ifTSl Eig. 4.1]. This function uses 
MATLAB’s built-in gmres function as well as the function zfmm2dpart from the fast 
multipole toolboxEMMLIB2D HU. The main inputs of the method consist of the discretized 
functions and 7 j(f) described above. 

Eor domains with smooth boundaries, we discretize the interval [0, 27r] by n nodes 
si,..., s„ and write s = [si,..., s„] where n is an even integer. For simplicity, we usu¬ 
ally take equidistant nodes, i.e.. 


Sk = {k- 1) — , 
n 

k = 1,... ,n. 

(3.6) 

Then ^ copies of s give a discretization 



t = [s, s,. . . , 

sf G 

(3.7) 

of the parameter domain J, leading to the discretizations 


?7(t) = [m{s),V 2 {s),. ■. ,r]e{s)f, vi^), 

7j(t) = - log|?7(t) - aji G C^”, 

(3.8) 


j = 1, 2,..., £. We store the discretized functions in the vectors et, etp, gam, respectively, 
and call 
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function mu = logcapacity(et,etp,alpha) 

% Computes the logarithmic capacity of a compact set with L components. 
% Input: 

% et = [eta_l; eta_L] (discretized boundary; column vector) 

% etp = derivative of the boundary curves (same format as et) 

% alpha = [alpha(1); alpha(L)] (auxiliary point alpha(j) interior 

% to j-th boundary curve Gamma_j) 

L = length(alpha); %% number of boundary components 
n = length(et)/L; %% number of nodes per boundary component 

% Auxiliary functions gamma_j(t) = -log |eta(t)-alpha_j| 
for k=l:L 

gamj(:,k) = -log(abs(et-alpha(k))); 

end 

% Compute the auxiliary functions h_j 

A = ones(size (et)); %% the function A in the gen. Neumann kernel 
for k=l:L 

[~,hjv(:,k)] = fbie (et, etp, A, gamj (:,k),n,5, [],le-14,100); 

end 


% Build and solve linear system for m_l, ..., m_L, 
for j=l:L 

for k=l:L 

hj(k,j) = sum(hjv(1+(k-1)*n:k*n,j))/n; 

end 


end 

matA 

matA(L+1,1:L) 
matA(l:L,L+l) 
vc_right 
vc_right(L+1) 


= hj; 

= 1 ; 

= - 1 ; 

= zeros(L+1,1); 

= 1 ; 


log(mu) 


X = matA\vc_right; 
mu = exp(x(L+1)); 
end 


Figure 3.1. MATLAB code for the computation of the logarithmic capacity. 


["',h] = fbie (et, etp, ones (size (et)), gam, n, 5,[], tol,maxit) 

Here [ ] means that GMRES is used without restart, toi is the convergence criterion used 
within GMRES, and max it is the maximal number of GMRES iterations. In the numerical 
experiments described in Section|4]we have used toi=le-14 and maxit = l 00. The output 
of fbie are the values hj{t) of hj from (13.41) , and the values hkj are computed by taking 
arithmetic means: 

^ kn 

^ n ^ tt): fc = l,2,j = l,2,...+. 

i—l + (fc—l)n 

These values are used to set up the linear algebraic system (13.51) . which we solve directly as 
mentioned above. 

Figure 13.11 shows our MATLAB implementation of the overall method described in this 
section, where the inputs et and etp are given as in (13.81 ). and alpha is the column vector 
of auxiliary points aj. 

The presented method can be extended to domains with corners. We assume that the 
corner points are not cusps and that the tangent vector of the boundary has only the first kind 
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discontinuity at these corner points. The left tangent vector at a corner point is considered as 
the tangent vector at this point. In this case, the solution of the integral equation ( 13.31 ) has a 
singularity in its first derivative in the vicinity of the corner points Using the equidistant 
nodes (13.71 ) to discretize the integrals in (13.31 ) and (13.41 ) yields only poor convergence lfT3l |22l 
flEi . To achieve a satisfactory accuracy, we discretize the integrals in (13.31 ) and (13.41) using a 
graded mesh which is based on substituting a new variable in such a way that the discontinuity 
of the derivatives of the solution of the integral equation at the corner points is removed. 

Following Kress ifTSll . we define a bijective, strictly monotonically increasing and in¬ 
finitely differentiable function, w : [0, 2tt] —>■ [0, 27r], by 


w{t) = 27r 


v{ty 


v{t)P + v{2tt — t)P ’ 


where 


v{t) = 



TT — t 


It — TT 


p TT 


1 

2' 


The grading parameter is the integer p > 2, and the cubic polynomial v is chosen to ensure 
that, for the equidistant mesh si, S 2 , • ■ •, Sn, almost n/2 of the grid points tu(si), w{s 2 ), ■ ■ ., 
w{sn) are equally distributed throughout [0, 27r], and the other half is accumulated towards 
the two endpoints 0 and 27r. 

Assume that the boundary T^ has qk > 0 corner points 


VkiO), rik{2Tr/qk), r]k{4:7r/qk), ..., rik{2{qk - l)Tr/qk). 


Then we define a bijective, strictly monotonically increasing and infinitely differentiable 
function, Sk ■ Jk ^ Jk, by 


Sk{t) = < 


w{qkt)/qk, 

{w{qkt - 2tt) + 2TT)lqk, 
{w{qkt - 47r) +47r)/gfc, 


t e [0,27r/gfc), 
t e [27r/gfc,47r/gfc), 
t G [47r/gfe,67r/gfc), 


yw{qkt - 2{qk - l)7r) + 2{qk - l)Tr)/qk, t G [2[qk - l)7r/gfc, 27r]. 


Since the function w has a zero of order p at the endpoints f = 0 and t = 2tt ifTSl Theo¬ 
rem 2.1], the function 5k is at least p times continuously differentiable. If the boundary T^ 
has no corner points we define the function 5k{t) by 


6k{t)=t, teJk- 


Finally, we define a function 5 : J ^ Jhy 


5i{t), t G Ji, 

b(i) = ; 

5e{t), tGJi. 

As noted in the first paragraph in lfT4l p. 242], using the graded mesh 


<5(t) = [<5i(s),J2(s),...,b,(s)]^GC^" 


(3.9) 
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for discretizing the integrals in ( 13.3b and ( 13.4b is equivalent to parameterizing the boundary F 
by r]{6{t)), and then solving the integral equation as in the case of smooth domains. Hence, 
we have the discretizations 

r]{6{t)), ri{5{t))S{t), 7 j(^(t)) = - log|?7((5(t)) - a^l e (3.10) 

j = 1.2.which replace the discretized functions in (13.8b . In our numerical experiments 

we have used the grading parameter p = 3. 

4. Numerical examples. We now present numerical examples that illustrate our meth¬ 
od. If not stated otherwise, we use an equidistant discretization of [0,27r]. Computations 
were performed in MATLAB R2013a on an ASUS Laptop with Intel Core i7-4720HQ CPU 
@ 2.60Ghz 2.59 Ghz and 16 GB RAM using the code shown in Figure [UT] Computation 
times were measured with the MATLAB tic toe command. 

4.1. Sets with one component. Example 4.1 (Disk). Let '■= {z & 'C ■. \z\ < r} 
denote the closed disk with radius r and c{Er) = r\ see Table l2.ll We use the parametrization 

p : [0, 27r] —>• dEr, t i-A re“'*, 

and n = 2® = 256 points in the discretization. For r = 1 our method computes the exact 
value c{Ei) = 1.0. For ?’ = 2 it computes the value 2.000000000000003, accurate to 14 
digits. (The relative error is 1.33 • 10“^®.) Each computation took less than 0.1 s. 

Example 4.2 (Ellipse). We consider the family of ellipses Ed with semi-axes a = 1 
and 6 = 10“^, and c{Ed) = (1 -f 10“‘*)/2; see Table 12. II A parametrization of the boundary 
is 


Vdit) = cos(t) — ilO ‘^sin(f), t £ [0, 27r]. 

We use d = 1, 2, 3, 4 and n = 2^ with fc = 8, 9,..., 18. For d = 1 we have c{Ei) = 0.55, 
and the value computed by our method has a relative error smaller than 10“^® for every n. 
For d = 2,3,4 the relative errors of the computed values are shown in Figure |4T| (left). We 
observe that our method is less accurate for larger d, i.e., for “thinner” ellipses. For large d the 
method still yields a very accurate approximation of c{Ed), but this requires a large increase 
of the number of discretization points. This may be related to the fact that the auxiliary point 
inside a “thinner” ellipse is necessarily closer to the boundary of the ellipse (cf. our comments 
after equation (13.2b ). On the right of Figure iTTI we show the computation times. 

Example 4.3 (Half-disk). Let E \= {z & \z\ < l,Im(z) > 0} be the upper half 

of the unit disk with c{E) = 4/3®/^; see Table l2.1] Since the boundary of E has corners, we 
use a graded mesh as described in Section[3 The following table shows the computed values 
(correct digits are underlined ), relative errors, and computation times for increasing n: 


n 

computed capacity 

relative error 

time (s) 

210 

0.769800347826294 

1.44- 

10"® 

0.1 

2“ 

0.769800357536609 

1.80- 

10-9 

0.2 

212 

0.769800358746887 

2.24- 

o 

7 

o 

1—I 

0.4 

213 

0.769800358897939 

2.80- 

10-11 

0.5 

214 

0.769800358916802 

3.51 • 

10-12 

0.8 

215 

0.769800358919097 

5.25 • 

10-1® 

1.5 

216 

0.769800358919667 

2.15- 

10-1® 

3.7 

217 

0.769800358919514 

1.63- 

10-14 

6.1 
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k with n = 2'' 


k with n = 2'' 


Figure 4.1. Results for Examole \4.2\ ( Ellipse): Relative errors of the computed logarithmic capacity (left) 
and computation times (in seconds, right). 


Example 4.4 (Square). Let E be the square with vertices 1 + i, 1 — i, —1 — i, —1 — i 
and hence side length h = 2, giving c{E) = ir(l/4)^/7r^'^^; see Table 12.1] Evaluating this 
expression in MATLAB gives c{E) = 1.180340599016096, which we use as the “exact” 
value for our experiment. As for the half-disk, our discretization uses a graded mesh. The 
following table shows the computed values (correct digits are underlined) , relative errors and 
computation times for increasing n: 


n 

computed capacity 

relative error 

time (s) 

2 ® 

1.180328330582103 

1.04- 

10-05 

0.1 

29 

1.180339089365394 

1.28- 

10-06 

0.1 

210 

1.180340411967884 

1.58- 

10-07 

0.1 

2“ 

1.180340575744745 

1.97- 

10-08 

0.2 

212 

1.180340596114299 

2.46- 

10-09 

0.2 

213 

1.180340598653831 

3.07- 

10-10 

0.5 

214 

1.180340598970863 

3.83- 

10-11 

0.6 

215 

1.180340599010508 

4.73- 

10-12 

1.1 

216 

1.180340599015215 

7.47- 

10-13 

2.3 

217 

1.180340599016100 

3.57- 

10-13 

5.9 


We also compute the logarithmic capacity with the Schwarz-Christoffel Toolbox ||7|. With 
the default settings the commands 

p = polygon ([ 1 + i,-1-l-i,-1-i, 1-i ]) ; 
f = extermap(p); 
capacity(f) 

yield the value 1.1803405990 90706, which has the relative error 6.32 • 10“^^. Eor a tolerance 
of 10“^^ instead of the default value 10“® the Schwarz-Christoffel Toolbox returns a very 
accurate result with the relative error 3.76 • 10“^®. 

In both Example 14. 3 1 and |4!^ the relative error in our method converges to the machine 
precision for increasing n. The reason that we need many points to obtain very accurate 
results is that in both examples the boundaries have corners. 

Example 4.5. We consider the set in Eigure 1472) (left), which is of the form introduced 
in na. Eor these sets an analytic parameterization of the boundary and the logarithmic 
capacity are known explicitly. Here we consider the compact set E bounded by 

77 (f) = fe[0,27r], 
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Figure 4.2. The sets from Examoles \4.5\ (lefi) and \4. SU risht). The (blue) dots show the auxiliary points ay. 


where tp is the conformal map given by HU Equation (3.2)] with the parameters Am = 1, 
(j) = 7r/2 and e = 0.1. Then c{E) = 1.223502096192244; see (El Corollary 3.4]. The 
following table shows the computed values (correct digits are underlined) , relative errors and 
computation times for increasing n: 


n 

computed capacity 

relative error 

time (s) 

28 

1.223385602611070 

9.52 • 10-°8 

0.1 

2^ 

1.223500703601890 

1.14-10-°® 

0.1 

210 

1.223502095786708 

3.31 • 10-1° 

0.1 

2“ 

1.223502096192245 

3.63-10-1® 

0.2 


As in the two previous examples, the relative error converges to machine precision as n 
increases. The reason that we need many points to obtain an accurate result is different 
from the two previous examples. Here the boundary is analytic, but with equispaced points 
in [0, 27r] and the above parametrization, only few discretization points lie on the inner arc, 
which is not well resolved for smaller n. Here, a different parametrization might lead to very 
accurate results already for small n, but we did not pursue this further. 

4.2. Sets with several components. Example 4.6 (Two disks with equal radii). Let 
r, zq G K with 0 < r < zq and let E be the union of the two disks Z)r(zo) = {z G C : 
\z — zo\ < r} and Dr{—zo) = {-2 G C : |z + zo| < r}. Erom (El Theorem 4.2] we know 
that 


where 


and 


ofc t - ._ 

c{E) = —Jzl - rV2i(l + C2), 

TT V 


\/zo + r — x/zq — r n TJ ( ^ ^ 

P= , 7^^^ L = 2pYY' 


y/zo + r + y/zo - r ’ 


fc=i 


8fc-4 1 ’ 


K = K{L^) = 


’o y/{l-t^)il-LH^) 


dt. 


The product for L converges very quickly, so that it suffices to compute the first few factors 
to obtain the correct value up to the machine precision. Using this value of L we evaluate the 
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complete elliptic integral of the first kind for K with the MATLAB command ellipk. We 
use the result as the “exact” value c{E). The constant L can also be written with the Jacobi 
theta functions as L = 02 ( 0 ;/ 5 ^)/ 03 (O; p^), using their product representation ll34l §21.3]. 
We apply our algorithm for Zq = 1 and r = 0.5, 0.7,0.9 with n = 2® (giving 2n = 2® nodes 
in total). The results are very accurate (14 digits of accuracy) and are shown in the following 
table (correct digits are underlined ): 


r 

computed capacity 

exact capacity 

relative error 

time (s) 

0.5 

1.030651235187015 

1.030651235187014 

8.62- 10-1® 

0.1 

0.7 

1.252472555601970 

1.252472555601971 

1.42- 10-1® 

0.1 

0.9 

1.465698640729795 

1.465698640729791 

2.42- 10-1® 

0.1 


Example 4.7 (Two disks with different radii). Let 0 < u < v and 


a = 


sinh(u) 
sinh(r; — u) 


and 


sinh(M) 
sinh(u — u) 


Then the capacity of Ea,r = Di{0) U Dr{a), which is the union of two disjoint disks, is 


c{Ea,r) = e“ sinh(M) 


02 (O;e-") 03 (O;e-") 04 (O;e- 


01 (m; e“") 


where 0i, 02 , 03,04 are the Jacobi theta functions (see below). This formula was brought to 
our attention by Thomas Ransford, who derived it from the Green’s function of an annulus 
in ||4] Ch. V]. We are not aware of a publication of this result in the literature. The Jacobi 
theta functions are ll^ §21.1] 


01 (z; q) = 2gi/4 sin((2n + l)z), 

n—0 
oo 

02 (z; q) = 2gi/4 ^ cos((2n + l)z), 

n—0 
oo 

(^ 3 {z',q) = 1 + 2 '^ cos(2nz), 

n—1 
00 

04 (z; g) = 1 + 2 ^(—cos(2 nz), 

n—1 

and we evaluate them by truncating the series when the absolute value of the terms become 
smaller than MATLAB’s eps. We use the computed value of c{Ea^r) as the exact capacity. 
We apply our numerical methods with n = 2® nodes per boundary, and obtain the very 
accurate results shown in the following table (correct digits are underlined) : 


{u,v) 

computed capacity 

exact capacity 

relative error 

time (s) 

(0.5,0.7) 

2.991271539541696 

2.991271539541695 

2.96- 10-1® 

0.2 

(0.5,1.0) 

1.637069166040759 

1.637069166040759 

2.72- 10-1® 

0.1 

(0.5,1.5) 

1.260209159232260 

1.260209159232259 

1.76- 10-1® 

0.1 


Example 4.8. Let E be the union of a disk and two half-disks, as shown in Eigure l4(2l 
(right). We are not aware of an analytic formula for c{E), but it has been shown in ll25ll that 


c{E) e [2.1969933, 2.2003506], 
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(b) 64 disks 
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Figure 4.3. The sets in Exanwle \4 . 91 The (blue) dots show the auxiliary points cxj. 


and the authors “best guess” is c{E) k, 2.19699371717. We parameterize the circle by 
771 (i) = and the two half-disks analogously to Example 14.31 With n = nodes per 
boundary component our method gives the computed value c{E) k. 2.196993710282112 in 
about 0.9 s. For n = 2^® our method computes the value c{E) 2.196993717171386 in 
36.2 s, and this value matches exactly the estimate from ll25l . 

Example 4.9. We consider the compact sets shown in Figure 1431 These sets were also 
used in 11211 Examples 2-5]. To our knowledge, the logarithmic capacities of these sets are 
not known analytically. The following table shows the values computed by our method and 
the computation times; 


E 

n 

computed capacity 

time (s) 

7 ellipses 

28 

4.961809958325545 

1.2 

64 disks 

28 

7.177814562549484 

43.6 

4 squares 

to 

0 

3.083190170261768 

1.3 

3 sets as in lfT2l 

to 

0 

2.977866214534663 

1.1 


Example 4.10 (“The World Islands”). We consider the unbounded domain K, of con¬ 
nectivity £ = 210 exterior to an artificial archipelago located in the waters of the Arabian 
Gulf, four kilometres off the coast of Dubai, and known as “The World Islands”; see also CD. 
An aerial image from Google Maps of “The World Islands” is shown in Figure llAl deft). The 
boundaries of the islands extracted from the aerial image are shown in Figure 14.41 (right). 
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(a) An aerial photograph of “The World Islands” (b) The boundaries of the islands extracted from the image 
Figure 4.4. “The World Islands” from Example \4.10\ The (blue) dots show the auxiliary points ay. 


and we parameterize them using trigonometric interpolating polynomials. The boundaries 
are very close to each other but they do not touch. The following table shows the values 
computed by our method and the computation times: 


n 

computed capacity 

time (s) 

2^ 

4.384226180107323 

484 

2^ 

4.388057916386704 

872 

2^ 

4.387882300144899 

1464 

28 

4.387881092385658 

2813 

29 

4.387881092740317 

4986 

210 

4.387881092740335 

9656 

2“ 

4.387881092740328 

17717 


4.3. Several real intervals and Cantor sets. In this section we consider sets consisting 
of several real intervals, and in particular the classical Cantor middle third set and generaliza¬ 
tions of it. Our method is not directly applicable to such sets, since these are not bounded by 
Jordan curves. We overcome this difficulty by considering a preliminary conformal map to 
“open up” the intervals and obtain a compact set of same logarithmic capacity, but bounded 
by smooth Jordan curves. 

Let E he a set that consists of i real intervals, so that its complement /C = = C\E 

is an unbounded domain (containing oo) bounded by i parallel straight slits, all on the real 
axis. In order to apply our method, we will hrst use conformal mappings to compute an 
unbounded multiply connected domain G exterior to a disjoint union of i ellipses, such that 
the domains G and /C are conformally equivalent. The computation of G and the conformal 
map z = uj{Q from G onto /C is based on a technique developed recently in which we 
describe in AppendixlAl The conformal map is normalized by uj(() = C + 0(l/() as ( —4 oo, 
which makes it unique, and also implies c(E) = c(G‘^) by ll23l Theorem 5.2.3], that is, the 
capacity of E can be computed as the capacity of the union of the £ ellipses. In summary, we 
use the method described in Appendix lAlfor computing G, and then apply our usual method 
(as stated in Fig. 13. lb in order to compute an approximation of c(E) = c(G‘^). 

We start with the much studied case of two real intervals, for which the logarithmic 
capacity is known analytically. The numerical results for these sets also give an indication of 
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the accuracy of our method for the Cantor sets, where no analytic formula for their logarithmic 
capacity is known. 

Example 4.11 (Two real intervals). Let Ea^b = [—1,«] U [b, 1] with —l<a<6<l. 
When the two intervals have the same length, i.e., when 0 < —a = b < 1, the exact capacity 
of Ea^b is given by (see Table l2.ll i 

c{Ea,b) = i\/l - a2. 

For the general case, an analytic formula derived by Achieser Cl (see also 1291) has the form 

(TP t — 1 +^ ^ 4 ( 0 ; g ) 

^ 2(i + a)04(AV2;g)’ 

where 

^ = k'=Vl-k\ K = K{k), K' = K{k'), q = 

K{k) is the complete elliptic integral of the first kind (see Example 14.61 above). 64 is the 
fourth Jacobi theta function (see Example l4.7l above). and 0 < A < AT is defined uniquely by 
the Jacobi elliptic function sn as 


sn(A, k) = 

We compute the values of K{k), K(k'), and 64 as explained in Examples 14.61 and 14.71 The 
value of the parameter A is computed using the MATLAB function asne, i.e., 

A = asne (x/(l-a)/2,fc). 

As explained above, we use our method to compute c(G'^) = c{Ea^b), where G is a 
domain bounded by smooth Jordan curves that is found by “opening up” the real intervals. 
Very accurate results are obtained with only n = 2® nodes per boundary, as shown in the 
following table (correct digits are underlined): 


(a, 6) 

computed capacity 

exact capacity 

relative error 

time (s) 

(-0.5,-0.1) 

0.488829271154718 

0.488829271154715 

4.77- 

10-1® 

2.9 

(0.5,0.6) 

0.499101557166360 

0.499101557166361 

1.11 • 

10-1® 

3.5 

(-0.5,0.3) 

0.457718411572721 

0.457718411572721 

8.49- 

10-16 

2.7 

(-0.5,0.5) 

0.433012701892217 

0.433012701892219 

5.64- 

10-1® 

2.1 

(-0.01,0.01) 

0.499974999374968 

0.499974999374969 

8.88- 

10-16 

3.6 


As a final remark concerning this example, it is worth mentioning that an analytic formula 
for the logarithmic capacity of sets consisting of several intervals has been derived recently 
in 13]. 

Example 4.12 (Cantor middle third set). In this example we consider the classical 
Cantor middle third set. Let Eq = [0,1] and recursively define 

Ek •= -^Ek-i U ^-Ek-i + 2 ^ > fc > 1. 

This means that Ek is constructed by “removing” the middle one third of each interval that 
Ek-i consists of. Then the Cantor middle third set is defined as E := While no 
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analytic formula for c{E) is known, several attempts have been made to numerically approx¬ 
imate c{E). In II 25 I it is shown that 

c{E) G [0.22094810685, 0.22095089228], 

and based on several of their computed values, Ransford and Rostand wrote that their “best 
guess” is 

c{E) « 0.220949102189507. 

Using two different approaches based on Schwarz-Christoffel mappings, Banjai, Embree, 
and Trefethen computed c{Ek) for k = 1,2,... ,9, and extrapolating from their computed 
values they obtained c{E) Ri 0.220949l|l Recently, Kruger and Simon lITSl obtained the 
value c{E) r; 0.22094998647421 in a study of the spectral theory of orthogonal polynomials 
associated to the Cantor measure. Referring to their result they noted that one “should only 
trust the first six digits or so”. 

We will now describe our approach for computing an approximation of c{E). Similar 
to Banjai, Embree and Trefethen, we will compute c{Ek) for a few small values of k, and 
then obtain an approximation of c{E) by extrapolation. We compute the capacities c{Ek) 
with the open-up method described in the beginning of Section l43] The resulting values and 
computation times for k = 1,2,...,12 are shown in the following table: 


k 

II 

to 

computed capacity 

time (s) 

1 

2 

0.235702260395518 

1.0 

2 

4 

0.228430704425426 

1.5 

3 

8 

0.224752818755436 

2.5 

4 

16 

0.222887290751916 

4.1 

5 

32 

0.221938129124324 

7.3 

6 

64 

0.221454205006181 

15.7 

7 

128 

0.221207178734289 

44.6 

8 

256 

0.221080995391656 

148.0 

9 

512 

0.221016516406108 

565.1 

10 

1024 

0.220983561713855 

2375.9 

11 

2048 

0.220966717159289 

9128.4 

12 

4096 

0.220958106742622 

34984.7 


In order to extrapolate from our computed values, we note that the differences 

dk = c(Ek) - c{Ek+i), k = l,2 ,...,ll 

behave linearly on a logarithmic scale; see the (blue) circles in Eigure l431 We therefore use 
the MATLAB command p=polyf it {l:ll,log{d{l:ll) , 1 ) ) for computing a linear 
interpolantp(a:) = pix -f p 2 of the values log(dfc). The computed coefficients are 

= -0.673356333942526, p 2 = -4.26116079806122, 

and the values p{k) for all fc = 1, 2,..., 48, where p(48) Ri 10“^®, are shown by the (black) 
pluses in Eigure l43] Since p{k) r; \og{dk), we can find an approximation of c{Ek) for each 
k = 13,14,... by extrapolation starting with our computed value for c{Ei 2 ), and obtain 

fc-i 

c{Ek) = c{Ei 2 ) - ^ exp(p(j)), k > 13. 

1=12 


^These computations, made in July 2005, were also reported in (25). 
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k 

Figure 4.5. The values = c{Ek) — for fc = 1, 2,..., 11 (circles) and the values exp(p(fc)) 

for k = 1, 2,..., 48 (pluses); see Examvle \4.I2\ 



Since p(49) < 10 we use only the values up to p(48), which gives our estimate for the 
capacity of the Cantor middle third set as 

c{E) K. 0.220949194629475. 

This estimate agrees up to the seventh digit with the estimates of Ransford and Rostand as 
well as Banjai, Embree and Trefethen. 

Example 4.13 (Generalized Cantor set). As in ll^ Section 6] we will now generalize 
the construction of the Cantor middle third set as follows. Let r G (0, 0.5) and Eq := [0,1]. 
Recursively define 


El := rEl_, U {rEl_, + (1 - r)) , k > 1, 

and let := The parameter r determines how much is removed from each interval 

of El_^ in order to obtain If we want to remove the middle q G (0,1), we need to set 
r = (1 — q)/2. Eor q = 1/3 we have r = 1/3 and hence E^^^ is the classical Cantor 
middle third set. The limiting cases are E^ = {0, 1 } with c{E^) = 0, and = [0, 1 ] with 
c(£;i/2) = 1 / 4 ; see Table O 

Using exactly the same approach as described in Example 14.121 we have computed the 
following approximations of c{E'^): 


q 

r 

computed capacity 

3/4 

1/8 

0.109156838696175 

2/3 

1/6 

0.13844418298159 

1/2 

1/4 

0.186511016338442 

1/3 

1/3 

0.220949194629475 

1/4 

3/8 

0.233218551525021 

1/5 

2/5 

0.23901897053678 

1/6 

5/12 

0.242233234580321 

1/7 

317 

0.244206003640726 

1/8 

7116 

0.245506481568117 

1/9 

4/9 

0.246410328817 

1/10 

9/20 

0.247064652445187 

1/11 

10/22 

0.247553947239903 

1/12 

11/24 

0.247929630663845 
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Figure 4.6. Computed approximation of c{E'") (circles) and the function f(r) (dashed); see Example \4. 1 3\ 


The (blue) circles in Figure l^hl show our computed approximations of c{E'^), where the 
right part of the figure is a closeup of the left part. The dashed line shows the function 

J.3 \ 3/2 

/(r) = r(l - r) - y - r j , 

which was suggested in as an approximation of c{E'^). The maximum distance between 
the values of f{r) and our computed approximations of c{E'^) is 7.5189 • 10“^. Thus, similar 
to computations reported in ll25l . the function /(r) very closely approximates our computed 
approximations of c{E^). 

5. Concluding remarks. We have presented a numerical method for the computation 
of the logarithmic capacity of compact sets bounded by Jordan curves in the complex plane. 
These sets may consist of several components and need not have any special symmetry prop¬ 
erties. In several numerical examples with sets for which the logarithmic capacity is known 
analytically, our method yields a computed approximation with a relative error close to the 
machine precision. For “simple” sets, in particular simply connected ones, the computations 
in MATLAB take at most a few seconds. 

Let us point out a few open questions. From a computational point of view, an automated 
choice of the auxiliary points aj (interior to each boundary curve) would be of interest. The 
numerical experiments for compact sets where the logarithmic capacity is known analytically 
suggest that the method is fast and accurate. A formal analysis of the numerical stability and 
accuracy of our method is beyond the scope of this article and remains a subject of further 
work. To compute an approximation of the capacity of the Cantor sets, we devised an ad hoc 
method to “open up” the intervals by conformal mapping and obtain a domain bounded by 
Jordan curves, to which our method could be applied. It would be of interest to devise an 
“open-up method” for general Jordan arcs that is computationally tractable. 

Appendix A. Numerical computation of the preimage of a parallel slit domain. 

Let n be a given parallel slit domain, i.e., the entire z-plane with m slits Lj, j = 
1,2,... ,i, along straight lines; see the top ofFigure lA.il An efficient numerical method for 
computing the conformal map 2 ; = u;(C) from an unbounded domain G exterior to i smooth 
Jordan curves Fj, j = 1,2,..., .f, onto the parallel slit domain fl such that w((^) = (^-l-0(l/(^) 
as 00 has been presented in ifTTl Section 4.5]. Assume that the boundary F of G is 
parametrized by the function r]{t) as in (13.11 1. Assume also the operators N and M are the 
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Figure A.l. The given parallel slit domain (top), the initial preimage domain (dashed line, bottom), 
and the computed preimage domain G (solid line, bottom). 


same operators as in (13.31) . Then we have the following theorem from ini. 

Theorem A.l. Let 

7 (f) = Im[77(f)], f e J, (A.l) 

let /i be the unique solution of the boundary integral equation 

(I - N)p = -M 7 , (A.2) 

and let h be the piecewise constant function 

h = (M/r - (I - N) 7 )/ 2 . (A.3) 

Then the function f with the boundary values 

fipif)) = 7 (f) + h{t) + i/r(f) (A.4) 

is analytic in G with f{oo) = 0 and the conformal mapping lo is given by 

c^(C) = C-i/(C), CeGur. (A. 5 ) 


In Theorem lA.il the domain G is assumed to be known, and the integral equation (IA.21 i is 
used to the find the conformal map z = uj(() from G onto the parallel slit domain ft = io(G). 
In our application with the Cantor sets, however, the domain G is unknown and the parallel 
slit domain 17 is known. Hence a straightforward application of a numerical method based on 
Theorem lA. 1 l is not possible. 

We will now describe an iterative method developed in ll20l for computing G and the con¬ 
formal map from G onto the (known) parallel slit domain 17. The method is an improvement 
of a numerical method suggested by Aoyama, Sakajo, and Tanaka 0, where the preimage 
G is assumed to be circular. Since the image region H is elongated (parallel slit domain), 
crowding can cause serious problems. Further, the convergence of the iterative method is 
slow if G is assumed to be circular. To overcome such difficulties, it was assumed in ll20ll that 
the boundaries of the domain G are ellipses instead of circles. 

Let \Lj I denote the length of the slit Lj and let Zj denote its center, j = 1,2,... ,i. In the 
iteration step i = 0,1, 2,... we assume that the domain G* is a multiply connected domain 
bounded by the ellipses T®, parametrized for j = 1, 2,..., £ by 

= Cj + 0.5(a® cost - i6® sinf), t £ Jj = [0, 27r]. 
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Then the following iteration computes the centers of the ellipses Q, the lengths of the major 
axes a*, and the lengths of the minor axes 6® for j = 1,2,. 

Initialization: 

Let e > 0 be a given tolerance and let Max he a maximum number of iterations. (In our 
numerical experiments in this paper we always used e = 10“^"^ and Max = 50.) Set 

Cj = Zj, a° = {1 - 0.5r)\Lj\, b° = ra°. 


where 0 < r < 1 is the ratio of the lengths of the major and minor axes of the ellipse (see 
Figure |AT| (dashed line, bottom) for r = 0.5). 

For i = 1,2,.. 

1. Map to a parallel slit domain 17® (based on Theorem lA.il) . which is the entire 
z—plane with i slits L®, j = 1,2,. .. ,i, along horizontal straight lines. 

2. If |L® I denotes the length of the slit L® and Zj denotes its center, then we dehne the 
parameters of the preimage domain G® as 


3. Stop the iteration if 


a® =a®-i-(|L®|-|L,|), 


(A.6) 

(A.7) 

(A.8) 


max (|z® - Zj\ + \\m - \Lj\\) < e or i> Max 

l<j<m 


Several numerical examples in this paper as well as in ||2l 123 show the convergence 
of this iterative method, but no proof of convergence has been given so far. Numerical ex¬ 
periments also show that the iterative method requires fewer iterations for small values of r, 
which means that the ellipses will be thin. For thin ellipses, however, we usually need a larger 
number of points n for discretizing the boundary integral equations and the GMRES method 
for solving these discretized equations requires more iterations to converge. 

In the numerical experiments with the Cantor sets shown in this paper we have not chosen 
to optimize upon these parameters, but we used the fixed values r = 0.5 and n = 64. The 
number of iterations for the convergence to the accuracy e = 10“^^ of the above iterative 
method applied in the computation of c{Ek) for k = 1, 2,..., 12 is shown in Figure lA!^ The 
(unpreconditioned) GMRES method for solving the discretized integral equations required 
between 5 and 11 iterations. 

Acknowledgments. We thank Nick Trefethen for sharing the numerical results on the 
capacity of the Cantor middle third set he obtained together with Banjai and Embree. We 
thank Thomas Ransford for bringing to our attention the analytic formula for the capacity of 
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